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•1  ABSTRACT 


The  steady  oscillatory  irrotational  motion  of  an  inviscid  incompressibl 
fluid  is  described  by  a  boundary-value  problem  of  elliptic  type.  A 
variational  form  of  this  problem  has  been  made  here  the  basis  of  a 
numerical  method.  The  problem  is  simplified  assuming  that  the  amplitudes 
of  the  generated  waves  are  small  compared  with  their  wave  lengths. 

The  numerical  satisfaction  of  the  radiation  boundary  condition  has 
been  investigated.  Some  sample  problems  with  known  solutions  have  been 
treated  first  in  order  to  test  the  method.  All  the  results  for  two- 
dimensional  motion  and  for  heaving  motion  of  an  axisymmetric  body  in 
infinite  or  finite  depths  show  very  good  agreement  with  existing  results. 
In  addition,  some  diffraction  problems  in  two  dimensions  with  homogeneous 
fluid  or  stratified  fluids  are  solved,  and  also  a  problem  with  non- 
uniform  depth. 

The  main  advantage  of  this  method  is  that  complex  geometry  of  the 
boundary  can  be  easily  accommodated;  for  example,  variable  depth  is  no 
more  difficult  than  constant  depth.  In  principle,  this  method  can  solve 
any  problem  of  elliptic  partial  differential  equations  with  boundary 
conditions  which  are,  partially  or  completely,  of  Dirichlet,  Neumann, 
or  mixed  type. 
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ABSTRACT 


The  steady  oscillatory  irrotational  motion  of  an  inviscid 

t 

incompressible  fluid  is  described  by  a  boundary-value  problem 
of  elliptic  type.  The  conventional  variational  form  of  this 
problem  has  been  made  here  the  basis  of  a  numerical  method  in 
which  an  appropriate  functional  is  minimized.  The  problem  is 
simplified  at  the  outset  by  replacing  the  nonlinear  boundary 
<  and  it  ions  by  linearized  ones  based  on  the  assumption  that  the 
amplitudes  of  the  generated  waves  are  small  compared  with  their 
wave  lengths. 

In  order  to  optimize  the  numerical  process,  the  decaying 
behavior  of  the  local  disturbance  has  been  investigated  and  the 
results  have  been  used  to  find  an  appropriate  position  for  imposing 
the  radiation  boundary  condition.  Some  sample  problems  with 
known  solutions  have  been  treated  first  in  order  to  cest  this 
new  method.  All  the  results  for  two-dimensional  motion  and  tor 
heaving  motion  of  an  axisymmetric  body  in  infinite  or  finite 
depths  show  very  good  agreement  with  extsting  results.  In 
addition,  some  diffraction  problems  in  two  dimensions  with 
homogeneous  fluid  or  stratified  fluids  are  solved,  and  also  a 
problem  with  non-uniform  depth. 

The  main  advantage  of  this  method  is  that  complex  geometry 
of  the  boundary  can  be  easily  accomodated;  for  example,  variable 
depth  is  no  more  difficult  than  constant  depth.  In  principle,  this 
method  can  solve  any  problem  of  elliptic  partial  differential  equa¬ 
tions  with  boundary  conditions  which  are,  partially  or  completely, 
of  Dirlchlet,  Neumann,  or  mixed  type. 
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Maximum  half  beam  of  cylinder. 

Submerged  area  of  the  cylinder  in  two 
dimensions. 

Water  plane  area  of  the  cylinder  of 
unit  length. 

Coefficient  matrix  and  its  components. 

Value  defined  in  (4-15). 

Body  boundaries. 

Complex  variables  c=a  +  ib,  and 
r  =  a  -  ib. 

Reflection  and  transmission  coefficients, 
respectively. 

% 

Decay  factor. 

\ 

Bottom  boundaries. 

General  functional,  cosine- 

and  sine-mode  functionals,  respectively. 

Functional  defined  in  an  element. 

Force  vector,  non-dimensional  force  vector 
and,  hydrostatic  force  vector,  respectively. 

Free-surface  boundaries. 

Free-surface  boundary  on  which  an  oscillatory 
pressure  distribution  is  specified. 

Acceleration  of  Gravity. 

Green's  function. 

Draft  of  a  floating  body  or  distance 
from  the  free  surface  to  the  top  of  a 
submerged  body. 

Defined  in  (4-14), 
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Depth  of  water. 
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Bessel  functions  of  tlie  first  kind,  of  the 
second  kind,  modified  Bessel  functions  of 
the  first  kind  and  of  the  second  kind,  all 
of  order  n,  respectively. 

Defined  in  (4-13). 

Eigenvalues. 

Moment  vector. 

Normal  vector. 

Interpolation  functions. 

Pressure,  its  cosine  mode,  its  sine  mode, 

and  the  nodal  values  of  the  cosine  and  sine  modes, 

respectively. 

Amplitude  of  the  pressure. 

Radiation  boundaries. 

Polar  coordinates  in  three  dimensions. 

Arc  length  on  the  body  profile. 

Time  variables. 

Period  (second). 

Velocity  components  in  the  x,  y  directions, 
respect ively. 

Velocity  vector,  its  cosine  mode  and  sine  mode, 
respectively. 

Cosine-  and  sine-moue  normal  velocities  and 
their  nodal  values,  respectively. 

Inertial  Cartesian  coordinates. 

Cartesian  coordinates  fixed  in  the  body. 

Amplitude  of  the  generated  waves. 

Wave  profile,  its  cosine-  and  sine-mode 
wave  profiles,  respectively. 

Non-dimensional  wave  profile,  its  cosine  and 
sine-mode  profiles,  respectively. 
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Diff racted-wave  profile,  its  cosine-  and  sine¬ 
mode  profiles,  respectively. 

Incoming-wave  profile,  its  cosine-  and  sine-mode 
profiles,  respectively. 

Amplitudes  of  the  incoming,  reflected,  and 
transmitted  waves,  respectively. 

Complex  variables  z  =  x  +  iy,  z  =  x  -  iy. 

Added-mass  and  damping  coefficients. 

Amplitude  of  sway  motion. 

Amplitude  of  heave  motion. 

Amplitude  of  roll  motion. 

Phase  angles,  for  the  diffracted  and  total 
waves,  respectively. 

Densities  of  the  fluids. 

Velocity  potential. 

Velocity  potentials  for  the  incoming  and 
diffracted  waves,  respectively. 

Velocity  potentials  corresponding  to  the  propagating 
wave  and  the  local  disturbance,  respectively. 

Cosine-  and  sine-mode  velocity  potentials, 
respectively. 

Local-disturbance  potential  due  to  a  source  in  ths 
fluid  and  the  image  of  the  source  in  the  upper 
plane,  respectively. 

Phase  lag  of  the  transmitted  wave. 

Tangential  vector. 

Circular  frequency  (radian/sec.). 

Infinitesimal  perturbation  parameter. 

Wave  number:  Vs- for  infinite  depth, 
and  **  V t&nk vH  in  water  of  finite  depth  H, 

Domain  of  definition. 
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Introduction 

When  we  treat  a  steady  oscillatory  irrotational  motion  of 
an  Inviscid  incompressible  fluid  described  by  a  boundary-value 
problem,  a  method  utilizing  Green's  function  is  used  most  often. 

But  sometimes  there  are  difficulties  in  carrying  through  compu¬ 
tations  which  are  involved  with  Green's  functions  and  in  addition 
the  method  has  another  drawback:  it  is  not  practical  for  a  very 
complicated  boundary  geometry,  i.e.,  for  example,  for  a  variable- 
depth  problem  ,  even  though  it  could  be  done  in  principle. 

In  this  paper  another  alternative  method,  a  variational 

method,  will  be  examined.  This  method  is  also  called  a  finite- 

< 

element  method  and  has  becom/e  a  very  useful  method  in  the  field 
of  structural  analysis  in  the  last  decade.  The  method  has  not 
been  used  much  in  other  engineering  fields, 

T^ere  do  exist  some  papers  which  treat  problems  of  fluid 
flow.  Zienkiewicz  (1964),  Matsuura  and  Kawakami  (1968), 

Zienkiewicz  and  Newton  (1969),  Holand  (1969)  and  Matsumoto  (1970) 
treat  oscillatory  motions  of  the  fluid  and  solve  for  the  pressure 
as  unknown  variable  in  a  finite  tank  since  the  pressure  is  the 
unknown  variable  in  their  formulations.  Taylor,  Patil  and 
Zienkiewicz  (1969)  treated  a  problem  of  undamped  harbor  oscilla¬ 
tion  based  on  the  shallow-water  theory.  Argyris,  Mareczek  and 
Scharpf  (1969)  and  Doctors  (1970)  treated  a  potential-flow  problem. 
Hunt  (1970,  1971)  treated  a  problem  of  sloshing  water  in  a 
container  based  on  a  disen  I  e-element  structural  theory  of  fluids. 

In  this  paper  the  behavior  of  a  local  disturbance  is 


t'xnm  inert  it:  considc rable  detail  and  the  radiation  condition 
is  examined  in  the  numerical  calculations.  The  minimizing 
functional  is  defined  as  a  function  of  the  velocity  potential, 
and  the  fir  si -order  linearized  problem  is  treated.  Several  two- 
dimensional  forced-motion  problems  and  three-dimensional  forced- 
heaving-mot ion  problems  of  an  axi-symmet rie  body  are  solved. 
Two-dimensional  diffraction  problems  are  also  solved  for  quite 
a  few  different  shapes  of  the  obstacles.  All  the  above  problems 
are  treated  in  water  of  finite  and  infinite  depth.  A  two 
dimensional  forced-motion  problem  is  also  treated  in  water  of 
variable  depth.  A  two-dimensional  diffraction  problem  in  the 
two  fluids  of  different  densities  is  treated. 

Most  results  are  compared  with  the  results  obtained  by 
the  other  methods  whenever  they  are  available.  Agreement  is 
generally  good. 

A  computer  program  has  been  written  that  can  solve  forced- 
motion  problems  or  diffraction  problems  in  a  homogeneous  fluid 
or  in  any  number  of  stratified  fluids  for  any  complicated 
boundary  geometries  in  two  dimensions.  Forced  heaving  motion 
of  any  axi-symmetric  body  in  three  dimensions  can  also  be  treated. 
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I.  MATHEMATICAL  FORMULATION 


We  assume  that  the  fluid  is  inviscid,  incompressible  and 
its  motion  irrotational ;  surface  tension  is  neglected.  Then 
there  exists  a  velocity  potential  which  governs  the  kinematics 
of  the  fluid.  Furthermore,  we  assume  that  the  motion  is 
sinusoidal  in  time,  so  that  we  can  drop  time  dependency  later. 
There  is  no  further  assumption  necessary  for  the  geometry  of 
the  boundaries. 

The  co-ordinate  system  is  right-handed  and  rectangular. 

The  y-axis  is  taken  directed  oppositely  to  the  force  of  gravity, 
the  x-axis  coincides  with  the  free  surface  when  the  fluid  is 
at  rest.  The  formulations  given  in  this  chapter  are  mainly  for 
the  two-dimensional  case,  but  one  can  find  the  formulations  for 
the  three-dimensional  case  in  the  reference  (Wehausen,  1960). 


1.  Governing  Equation  and  Boundary  Conditions 

When  we  define  a  velocity  potential  i£(x,y,t)  in  the 
fluid,  this  satisfies  a  Laplace  equation  throughout  the  xluid: 


V*  1  C  =  0  , 


(1-1) 


where 


> 


9 


VU.jj.fc) 


(1-2) 


In  order  to  complete  the  free-surface  boundary  conditions. 


it  is  necessary  t  *  make  use  of  one  >ther  < condition  beside^  *‘he 
kinematic  condition  on  the  free  surface,  since  the  pressure  is 
proscribed  but  the  form  of  the  surface  is  not  prescribed  a  priori. 
Thus  wo  use  both  kinematic  and  dynamic  boundary  conditions  on 
the  free  surface.  From  them  we  obtain  on  the  free  surface, 
y  -  Y(x,t),  the  boundary  conditior 

+  +  +5,^*0  .(1-3) 


For  a  moving  boundary  the  kinematic  boundary  condition  is 

§f-  =  n-v< ,  d-4) 

where  tl  is  a  normal  vector  to  the  surface  and  V(x,y,t)  may  (>e 
taken  as  known.  When  the  boundary  is  fixed  in  space,  the  equation 
(1-4)  degenerates  into  the  homogeneous  boundary  condition 


jtM-cx,  y,fc) 

a-n. 


=  O 


(1-5) 


When  the  depth  of  the  fluid  is  infinite,  we  have 


jLurvx- 

JT-*0 


(1-6) 


'"he  radiation  condition  requires  the  waves  to  be  progressing 
outwards  from  the  wave-generating  source  and  imposes  a  uniqueness* 


*  A  partial  differential  equation  of  elliptic  type  with  boundary 
conditions  of  Dirichlet,  Neumann  or  mixed  type  requires  the  boundary 
to  be  closed  to  ensure  a  unique  and  stable  solution.  See  Morse  and 
Feshbach  (1953),  pp,  70(5. 
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which  would  not  otherwise  be  present. 

Referring  to  Fig.  1,  the  cylinder  is  described  parametrically 
by  Cartesian  coordinates  x  =  x(s),  y  =  7(s)  in  a  system  fixed  in 
the  cylinder  such  that  this  coordinate  system  coincides  with  the 
inertial  coordinate  system  Qxy  only  when  the  body  is  at  rest. 

The  forced  oscillation  of  the  body  is  defined  by  giving  the 
coordinate?  of  the  origin  of  the  body  coordinate  system  and  the 
angle  betwoen  Ox  and  57: 

|(t)  =s  L$,ct>  »£.  C'l,  COS  ft  t  U  sJ—r t)  J 

» t  %<t)  c*»  ft?  t  flwrt)  t  (1-7) 

eib)  a  £  ftfc)  -  tc  o,  u*  +  A  frn.fb)  , 

where  7<-  7*.  A.  0i  are  constants  and  $£fc>  ,  ,  0fc) 

describe  sway,  heave,  and  roll  motion,  respectively.  The 
infinitesimal  perturbation  parameter  g  measures  the  relative 
size  of  body  motion  and  hence  of  the  motion  of  the  fluid 
throughout  the  domain. 


Fig.  1  Coordinate  Systems 
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(1-8) 


Th('  motion  of  a  point  of  (he  surface  of  the  body  described 
in  terms  of  ihe  inertial  coordinates  is  as  follows: 

XCs.t;  *  +  t  [  -  ?«>  &*)]  +  Oli1)  , 
y  CS.  fc)  =  ys>  +  l  [  +  XU  fa « J  +  0&‘>  , 

and  the  normal  velocity  of  the  body  surface  is  given  by 

»•*  «  H *&*9  *  00  (1-.) 


on  (  x(s,t),  y(s,l)). 

The  boundary  condition  (1-9)  is  not  yet  a  tractable  boundary 
* 

condition  ,  unless  one  solves  this  problem  at  the  instantaneous 
position  or  as  an  initial-value  problem  because  the  boundary 
geometry  is  not  fixed  in  space  but  moving  as  a  function  of  time 
t.  Hence  the  domain  of  definition  in  which  a  Laplace  equation 
is  to  be  satisfied  is  not  fixed  but  is  also  moving  as  a  function 
of  time.  To  solve  this  as  an  initial -value  problem  is  out  of  this 
paper's  scope,  therefore  one  should  try  to  obtain  a  tractable 
boundary  condition  for  a  boundary  fixed  in  space  for  all  time. 

In  order  to  express  a  normal  velocity  distribution  on  a  boundary 
fixed  in  space,  we  shall  make  use  of  the  analytic  continuation  of 
the  velocity  potential  <£(x,y,t),  so  that  the  potential  evaluated 
on  the  body  can  be  expanded  in  Taylor  series  in  £,  at  a  reference 


*  It  is  of  interest  to  note  that  the  differential  operator  for  a 
Laplace  equation  does  not  have  an  operation  with  respect  to  time 
whereas  a  moving  boundary  is  a  function  of  time.  A  boundary  con¬ 
dition  on  a  moving  'material'  boundary  is  always  described  in  a 
Lagrangean  representation  rather  than  an  Eulerian  representation. 
However,  in  a  diffraction  problem  with  a  fixed  body  in  incoming 
waves,  the  boundary  condition  on  the  body  is  Eulerian.  (See 
Kq.  (1-25). 
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boundary  (usually  the  average  position  Is  taken)  fixed  in  space 
for  all  time  and  then  in  a  perturbation  series  in  £,  in  terms  of 
the  field  variables: 


*  0(i)  ,  <1“10> 

3E  <*.*.+>  +6* 2**fe a (i*’*i> 

We  obtain  (1-12)  for  the  normal  derivative  of  the  potential 
evaluated  on  the  surface  Sq,  i.e.,  (  x=x(s),  y=y(s)  )  :  N 

~  ik  $'k.3.v  +  cxti  <i-») 

Equations  (1-9)  and  (1-12)  yield  the  required  boundary 
condition  on  the  body  (1-4)  for  a  forced  motion. 

Next  we  express  the  free  surface  Y(x,t)  in  a  power 
series  of 

*  &  Y  tx.h)  -t  E,lY^lz't}+  .  .  •  (1-13) 


Following  the  same  procedure  used  to  transfer  the  boundary  condition 
for  a  moving  body  to  a  fixed  boundary,  we  shall  expand  j£(x,Y,t) 
in  a  Taylor  series  about  the  neutral  free  surface  y  =  0  and 
substitute  (1-13)  in  (1-3).  Then  we  will  obtain  the  following 
first-order  free-surface  boundary  condition: 


—  ■x?ul  (1-14) 

*,*.<->  *  tZ,  t  u  l4) 

When  we  drop  the  superscript  in  the  potential  for  convenience 

and  rewrite  as  a  first-order  problem,  we  have 

V1  3?C*,d.fc  >  *  0  (in  the  fluid) 


<£tv  +  9-  -  o 


(on  y=0) 


5-n.  *  " 


(on  S  ) 
o 


(1-15) 


Further,  we  assume  that  the  potential  ^(x,y,t)  can  be 


decomposed  into  n  cosin. '-mode  potenti  i  and  a  sine-mode 
potential  as  follows: 

=  fVf)  U> s  &  +  f  sr-n,  Cb  #  (1-16) 

From  (1-7)  and  (1-9),  we  obtain 

v*  =*• 

(1-17) 

V*  (x,!).*)  *  £lC\x+%6i)  coifh  -(%  -.  Zb)  ffartjj 

where 

V  =  (V,  ,  tf,  ) 

From  (1-15),  (1-16)  and  (1-17)  we  obtain 


x7l  (ft(z.y>*o  , 

(1-18) 

t,  -  f  f  • 

f  f  (on  2*°)  * 

(1-19) 

(err.  *.  / b»dj )  t 

=* 

(  (rr.  •=-  ~W.  /*r  H) 

(1-20) 

^.-•o  3  > 

(1-21) 

8S— v  . 

f£  -»./ 

s*n.  > 

(1-22) 

where 

V(*.^,fcl  «  005  Ft  b  Vra,^>  sr-n.  Tt  J 

=  ^  -*  A?,  ro?t*-z&) )  , 

(1-23) 

vf=*  ( <rr-*,*y^, -re?,  *■**,» , 

As  the  radiation  condition  (Wehausen,  p.  480;  Stoker,  p,  174) 


we  can  assume  the  asymptotic  behavior: 


£  v^s  “*  * 
?X  £  v  f C  — »  0 


0A>  X 


> 
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X  oo 


(1-24) 


where  V  is  the  wave  number^i.e.,  for  infinite  depth  and 

■*  y  ficnJi  VH  for  finite  depth  H. 

When  we  want  to  solve  a  diffraction  problem  due  to  a  fixed 
body  in  the  fluid  with  a  homogeneous  incoming  wave,  we  assume 
that  the  total  potential  can  be  decomposed  into  a  known  potential 
of  an  incoming  wave  JjS^Xjyjt)  and  a  potential  of  a  diffracted 
wave  <§0(x,y,t).  Then  the  boundary  condition  on  the  body  for  the 
di f f racted-wave  potential  is 


fxOc.y.+)  . 


(1-25) 


Let  us  next  suppose  that  we  have  two  fluids  of  densities 
f,  ,  fi  ,  one  beneath  the  other,  the  common  surface  (when 
undisturbed)  being  plane  and  horizontal,  as  shown  in  Fig.  2, 


When  we  define  the  velocity  potentials  5,  .  &  in  the  upper 
and  lower  fluids,  respectively,  and  ".hen  we  use  both  kinematic 
and  dynamic  boundary  conditions  at  the  interface,  we  can  write 
the  first-order  linearized  boundary  condition  at  the  interface 


of  the  two  fluids,  i,e.,  y  =  -h,  as  follows: 

e,  *,«-?,  I* +&<>%*.,  = o  . 

e,  !>«.  +(e,-C)  j  =  <? . 


When  we  define 


3>.  «  £c  sr-n.  rtr  ,  t»>,  a., 


(1-26)  becomes 


1 

(&<)%■< 

crk-e,)j 

*  «&-  J 

?c  —  —  f  ^  -  <A 

*•  " 

+£JL- 

y5-  £lS— 

<6  <) » 

* 

(1-26) 


(1-27) 


(1-28) 


(1-29) 


It  should  be  noted  that  (1-28)  and  (1-29)  are  nuxed-type 
boundary  conditions  in  the  upper  fluid  and  the  lower  fluid, 
respectively,  and  that  thty  are  coupled  to  each  other. 
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When  a  pressure  distribution  is  confined  to  a  segment 
of  the  free  surface,  i.e., 


fCX.tr)  = 


f s  sr-rv  <r-b  ,  1*1  . 

o  ,  |*|  yl  . 


(1-30) 


with  y*  a  constant,  we  obtain  the  following  boundary  condition 
on  the  free  surface  from  (1-16)  and  (1-30): 


K  =  ° . 


1*1  **  , 


2*0, 


1*1  1A.  , 

-oo  <  X  <  -too  t  tf-Oj 


(1-31) 


where 


6  *  -  Ps*/e} 


Three  dimensions  will  be  considered  very  briefly  here. 

When  we  consider  the  three-dimensional  problem  analogous  to 
(1-18)  -  (1-31)  for  two  dimensions,  we  can  still  use  most  of 
the  two-dimensional  formulations.  We  take  the  velocity  potential 
$CX.  $,*,■&)  in  rectangular  coordinates  or  $  C  in 

cylindrical  coordinates.  The  Laplace  operator  becomes 


or 


where 


j.  3  /  a  A—  -1—  i—  -f-  — - 

V  -  T  ***■'  **■?#<*  * 


K  3* 

X  R  cos  <*, 


3-  R  SP*.«4 

In  (2-22)  we  now  have  VC= VCCX.<j.g) ;  *  and  0- -23)  will 

take  a  different  form  because  the  given  motion  has  now  six  degrees 


of  freedom. 
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As  the  radiation  condition  in  three  dimensions,  we  assume 
the  following  asymptotic  behavior: 


jr  <If  -;H *■>•<>, 


(1-32) 


where  , 

=  <fc  005  +  jP*  w-n  rfr, 

-  fVw  *  *■  f  sc*w 


or  equivalently 

ifw  +vrs>  -o  , 

im  -vfcJ  =0. 


(1-33) 


(1-34) 
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2.  Pressure.  Force.  Moment,  and  Wave  Profile 


The  velocity  potential  J(x,y,t)  defined  earlier  does 
not  give  any  dynamic  properties  of  the  flow  unless  Euler's 
Integral  Is  utilized.  The  pressure  jo  is  determined  by  Euler's 
integral , 


(1-35) 


with  the  term  suppressed  as  being  of  higher  order. 

The  first  term  in  (1-35)  is  a  hydrodynamic  pressure  and  the 
second  a  hydrostatic  pressure.  Then  the  force  on  a  body  submerged 
in  the  fluid  is  given  by 


p  ~  j  r « 
& 


(1-36) 


and  the  moment  with  respect  to  the  origin  by 


M  =  [  Pl'Xtlds  , 
f>o 


(1-37) 


where  ,S#  describes  the  interface  of  the  body  and  the  fluid  and 
tl  is  the  normal  vector  into  the  body. 

If  the  free  surface  is  described  by  y  =  Y(x,t),  then 


Y  “  to 


(1-38) 


When  we  decompose  Y  into  a  cosine  and  sine  mode  in  time. 


Y<7.-e)  =  Y<*>  C*S  <Tt-  +  YS<»  5fVi  6’-tr 


(1-39) 


thon  we  obtain,  from  (1-16),  (-36)  and  (1-39) 


Yw=  -f 


(1-40) 


When  a  body  is  undergoing  a  forced  oscillation,  it  is  more 
convenient  to  use  non-dimensional  quantities  for  (1-35)  - 
(1-40).  If  we  assume  the  forced  motion  has  only  heave  motion 
with  amplitude  V  ,  then  we  can  define 


?  -  */«?  , 

f  - 

7  =  y./7  ,  T-  Y'/7  ,  Y‘  -  YVf  , 


(1-41) 


where  A  is  water  plane  area  and  ,  F  ,  Y  are  the  amplitudes 
w  o  o’  o 

of  the  pressure  1*  ,  of  the  force  p  and  of  the  generated  waves  Y, 


respectively,  these  being  sinusoidal  in  time  -fc 


The  added  mass 


h 


*  * 
and  damping  coefficient  can  be  defined  as  follows  : 


a 

A 


(1-42) 


n  —  (n  i ,  n»> 


*  The  notation  JXai  and  has  been  introduced  here  In  order  to 
conform  with  the  more  general  notation  used  when  all 

degrees  of  freedom  are  present.  Such  problems  can  also  be  handled 
by  the  method  to  be  described. 
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where  A  is  the  subjerged  area  of  the  body  in  two  dimensions. 

One  can  define  all  these  non-dimensional  quantities  and  hydro 
dynamic  coefficients  for  other  modes  of  motion  in  two  dimensions 
or  in  three  dimensions,  but  these  will  not  be  given  here. 
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II.  LOCAL  DISTURBANCE,  PROPAGATING  WAVE  AND  RADIATION  CONDITION 


In  this  chnpter  an  attempt  is  made  to  obtain  a  more  tractable 
form  for  the  radiation  condition  in  order  to  facilitate  applying 
a  numerical  method.  One  main  difficulty  in  the  radiation  condi¬ 
tion  is  that  it  should  be  applied  as  x~*oo  (  or  R -*«  ).  But  one 
;  cannot,  of  course,  go  to  infinity  and  an  infinite  boundary  has 

I 

I 

always  to  be  truncated  at  some  ’  sufficiently  large'  distance  in 

J 

I  a  numerical  method.  In  this  chapter  some  criteria  are  developed 

» 

to  determine  how  near  to  the  moving  body  (or  a  fixed  body  in  a 
,  diffraction  problem)  one  can  construct  an  imaginary  boundary  on 

which  one  can  apply  a  more  tractable  version  of  the  radiation 

;  ‘  condition.  The  two-dimensional  case  will  be  treated  first, 

( 

i  In  order  to  examine  the  behavior  of  a  local  disturbance,  we 

! 

make  use  of  two  fundamental  tools  in  potential  theory.  The  first 
is  the  eigenfunction  expansion,  and  the  second  is  the  potential 
of  a  pulsating  source,  since  the  source  is  the  slowest-decaying 
singularity  among  all  orders  of  singularities.  It  seems  to  be 
easier  to  use  a  pulsating  source  in  the  case  of  infinite  depth, 
whereas  both  can  be  used  in  the  finite-depth  problem  and  give 
the  same  result. 

In  two-  and  three-dimensional  cases  much  is  very  similar. 
Therefore  in  section  two  of  this  chapter,  treating  the  three- 
dimensional  case,  details  will  be  mostly  omitted. 
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1.  Two-dimensional  Problem 

a)  A  Local  Disturbance 

In  order  to  make  use  of  the  eigenfunction  expansion  in  the 
finite-depth  problem,  we  construct  two  imaginary  vertical 
boundaries*  which  e> tend  from  the  bottom  to  the  free 

1.  Ct 

surface,  Including  any  moving  body  in  between.  First  we  assume 
the  depth  is  constant  at  y  t  -H.  One  can  subdivide  the  region 
52.  bounded  by  y=0,  y=-H,  Si,  and  %%  into  three  subregions 
Sit  ,  Slx  tSl^  as  shown  in  Fig.  3  so  that 

SI  =  SI,  *Si»  .  (2-i) 

Once  we  have  solved  the  original  problem  with  given  motion 
on  the  moving  body,  we  know  the  solution  5<*,y  ,t)  everywhere  in 
the  fluid.  Suppose  the  function  ?  «£(x,y ,t)/ >n  on  >8>2  is  computed 
from  this  solution  5.(x,y,t).  With  this  we  may  now  solve  a  new 
boundary-value  problem  in  the  subregion  52*  with  the  boundary 
condition  on  \£>^  derived  from  the  original  formulation.  We  should 
obtain  a  solution  identical  with  the  original  solution  in  the 


*  One  can  also  construct  only  one  vertical  imaginary  boundary 
which  pierces  the  body  and  extends  the  original  domain  of 
definition  to  include  the  immersed  part  of  the  body.  Then  we 
separate  the  original  problem  into  two  problems  for  which  the 
domains  extend  to  infinity  on  the  left-hand  side  and  right-hand 
side,  respectively.  One  assumes  that  there  exist  normal- 
velocity  distributions  for  each  sub-region  which  represent  the 
body  motion  with  ,ri  appropriate  restriction  on  the  body  geometry. 
In  the  three-dimensional  case  one  can  construct  an  imaginary 
vertical  circular  cylinder  which  contains  the  body  and  extends 
to  the  bottom. 
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Fig.  3.  Imaginary  Boundaries,  ifi/  and 

sub- region  Six  due  to  the  uniqueness  theorem.  A  similar 
procedure  can  be  applied  for  the  sub-region  SI)  . 

Let  us  consider  solving  the  problem  in  the  sub-region 
where  the  boundaries  consist  of  a  free  surface  $2*  the 
tious  moving  boundary  $  ,  the  bottom^,  and  the  boundary 
(x  =  oo  ).  We  may  shift  the  y-axis  without  loss  of 
generality,  as  shown  in  Fig.  4. 


Fig.  4.  Normal  Velocity  Distributions 
on  an  Imaginary  Boundary 


Now  this  new  problem  can  oe  interpreted  as  a  flexible-wall 
wave-maker.  In  order  to  solve  this  problem  one  can  use  a 
classical  method,  i.e.  the  separation  of  variables.  By  this 


method  one  can  easily  obtain  the  eigenfunctions  (Wehauaan, 
pp.  472-475).  The  eigenfunctions  are 
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[  •yti’oC.y+H),  oos  ini  #  (2-2) 

whe  re  , 

6*Vj  -  tw,  -totm-k 

..  .  (2-3) 


Fig.  5  shows  the  relations  (2-3)  for  the  eigenvalues,  and  the 
functions  (2-2)  may  be  shown  easily  by  direct  computation  to  be 
orthogonal  on  the  interval  -|j  S  ^  S  o ,  Both  orthogonality  and 
completeness  are  consequences  of  the  Sturm-Llouville  theory. 

By  making  use  of  these  eigenfunctions,  one  can  obtain  the 
solution  of  this  problem  as  follows: 

<£cx.^/b)  as  OOi  **  +  f*  €4: 


&SL  C*sh  Si'-n.  (tHoX  -  <2_4> 


o*  a}  -mix 

-Z.  -=rC  CCS  Y*ity+ll)  eaifltr 


where 


sr^.k  2m,H Ja  '> 

/’*»  m* <2-5> 


c  t  » 1,  ■  ) 
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Fig.  5.  Eigenvalues 


Further,  the  potential  3L  can  be  decomposed  Into  two  parts: 

one  represents  a  propagating  wave  system,  the  other  a  pure  local 
* 

disturbance  ,  i.e. , 

<g  CX.  -  Jp  +  t,  «  *.*)  ^  (2-6) 

where 

Spd.y.V  =  WWtJ-t-rt?  **•  (**-«*),  (2_7) 

c$Ct.  (2-8) 

Let  us  consider  the  behavior  of  the  local  disturbance  expressed 

in  the  equation  (2-8).  How  far  should  one  go  along  the  x-axis 

in  order  to  have  the  local  disturbance  reduced  to,  say,  0.5%  of 

the  maximum  (x=0)?  To  determine  this  we  have  to  examine  the 

*  One  can  also  obtain  the  solution  for  the  local  disturbance 
£.  (x,y,t)  by  imposing  £(x,y,t)-*.  O  ns  x-»*»  .  In  this 

case  one  will  obtain  {  c*s  ty-Hoj-  ,  i=l »2, ’ * ’  ,  as 

eigenfunctions. 


values  of 
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Mi,  (1=1,2,  •••).  The  value  m^H  is  the  smallest 
atyong  m^H  (1=1,2,* ••)  since  is  the  smallest  eigenvalue. 

From  Fig.  5  the  following  inequalities  hold: 

<  W|  H  <  } 

3 yx  <mxt\  < 

:  (2-9) 


Let  us  define  the  decay  factor  from  (2-8)  by 


d  5 


-WiX 


(2-10) 


When  we  non-dimens ionalize  the  variable  x  by  taking  K  as  a 
non-dimensionalizing  length  scale,  we  have 


<1  = 


-WiHX 

<? 


(2-11) 


where 

X  -  X/H 


(2-12) 


From  (2-9)  and  (2-11)  one  obtains 


e 


Kx 


i  <  e 


(2-13) 


Now  one  can  see  that  for  finite  depth  the  decay  factor  is 
independent  of  the  frequency  of  the  motion,  i.e.  of  &  or  *3C%!k\ 
For  example,  when  one  wants  to  assure  J  <  e"***  ,  then  one 

must  have  %  ,  that  is,  x  should  be  at  least  four 

times  the  depth.  When  one  considers  a  particular  frequency  such 
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that  the  wave  length  is  four  times  the  depth  H,  then  the  local 

disturbance  decreases  to  S  times  its  original  size  in 

a  half  wave  length  away  from  the  wave  maker.  Another  extreme 

example  is  the  case  in  which  the  frequency  of  the  motion  is  so 

high  that  the  wave  length  is  a  tenth  of  the  depth  H.  Then  the 

— nc> 

local  disturbance  decreases  to  g  times  its  original 

size  in  forty  wave  lengths  away  from  the  wave  maker.  Fig.  6 
shows  qualitatively  these  two  extreme  cases. 


It  was  mentioned  at  the  beginning  of  this  chapter  that  an 
estimate  of  the  extent  of  the  local  disturbance  could  also  be 
obtained  from  a  source  of  pulsating  strength.  For  the  case  of 
infinite  depth  we  shall  use  this  alternative.  A  source  with  pul¬ 
sating  strength  cos  at  (a.b)  is  given  (Wehausen,  p.  481) 

as  follows: 


S~t t. 


'IV  CZ-t) 


-  e 


(2-14) 
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where 


2  =» 

X  ^  , 

C  * 

CL  + , 

C  = 

fl-  -b  t  , 

V  = 

After  making  use  of  a  contour  integral  in  the  complex  plane 
and  some  manipulations,  one  can  reduce  (2-14)  to  (2-15)  and 
obtain  (2-16)  from  (2-6),  The  procedure  for  deriving  (2-15) 
from  (2-14)  is  given  in  Appendix  A, 


$p(xJ  a. vCx-*)  -ft] 


(2-15) 


where 


2  -<L 


2-C 


e.e 


u9i 


Now  one  can  examine  a  local  disturbance  defined  in  the  first 
equation  of  (2-15)  and  one  can  utilize  the  Gauss-Laguerre  formulas 
[Abramowitz  and  Stegun,  1967]  for  the  numerical  integrations. 

One  can  further  decompose 


+  , 


where 


-  f- 


-  -fv(i 

-  sr-K.  $>»;  e 
tg.- «r-n.  *  £»ift 


(2-16) 


(2-17) 


(2-18) 
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; 

J>LS  is  the  potential  due  to  the  source  and  is  due  to  the 

image  of  the  source,  i.e.  a  sink  in  upper  half-plane.  Xn  the  case 
of  infinite  depth,  one  can  obviously  define  the  decay  factor,  but 
'  Figures  7  to  10  show  the  values  of  the  potential  for  a  local 

disturbance  directly;  this  gives  enough  information  about  its 
behavior. 
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CE  DUE  TO  AN  IMAGE  CP 


THE  LOCAL  DISTURBANCE  DUE  TO 
PULSATING  SOURCE  AT  Y=-10V  C  Total)- 
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b)  Propagating  Waves  and  Habitation  Condition 

From  the  previous  sub-section,  one  sees  that  the  velocity 
potential  can  be  represented  as  a  sum  of  a  potential  of  a  local 
disturbance  and  a  potential  of  propagating  waves.  This  decomposi¬ 
tion  is  unique.  Now  we  can  examine  a  pure  propagating  wave  which 
is  the  solution  for  a  homogeneous  boundary  condition  with  no 
body  in  the  fluid  and  with  uniform  constant  depth  or  infinite 
depth.  The  solutions  are  given  in  many  books  on  hydrodynamics 
(e.g. ,  Lamb  p,  369).  They  are  as  follows:  .  • 

for  finite  depth, 

£  CtCj  -  A  ^>sk  +  toltyiL-Gt)  ,  (2-19) 

where 

^  a=  wave  amplitude 

for  infinite  depth 

35^  -*  3  eV*  c <05  (yx-ffc)  J  (2-20) 

where 

ys  f  Q.  *s  wove  amplitude 

In  the  formulation  in  Chapter  One,  the  radiation  condition 
simply  states  the  asymptotic  behavior  of  the  generated  waves  at 
a  large  distance  from  the  moving  body,  but  for  a  numerical  method, 
we  need  a  more  tractable  boundary  condition  than  (1-24)  or  (1-34), 
which  simply  state  an  asymptotic  behavior  at  infinity.  Therefore 
in  this  subsection  we  will  try  to  obtain  a  more  tractable  radiation 
condition  from  the  examination  of  pure- propagating  waves  in  order 
to  apply  it  to  the  numerical  scheme. 
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Let  us  try  to  derive  a  boundary  condition  from  (2-19)  or 
(2-20),  First  (2-19)  and  (2-20)  are  slightly  changed  in  (2-21) 
bv  defining  which  correspond  to  the  cosine  and  sine 

modes,  respectively: 

^  —  <j> *  coS-G*  +  $T-n.€t  (2-21) 


where 

for  finite  depth, 

yC  -  A  , 

Vs  -  A  Cssh  SPnVZ  , 

for  infinite  depth, 


fl 

=  b 

«** 

L+S  vx 

ioS 

•es  £ 

v? 

e 

S~-y^  VX 

(2-22) 


(2-23) 


Since  the  formulation  in  Chapter  One  is  given  for  the 
potentials  y’fx.jjland  after  the  time  has  been  precipitated 

out,  one  obviously  may  try  to  make  a  condition  (a  coupled  relation 
between  Jf**  and  )  from  (2-22)  or  (2-23),  One  readily  sees 

that  a  coupled  relation  cannot  be  made  from  them  unless  the 
derivative  terms  are  taken  into  consideration.  The  first  derivatives 
with  respect  to  x  and  y  are  as  follows: 
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for  finite  depth, 
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€  57-in.  VTL 

• 

(2-24) 


(2-25) 


From  (2-22)  and  (2-24)  or  (2-23)  and  (2-25)  the  following  coupled 
first-order  differential  equations  can  be  obtained  for  finite  or 
infinite  depths: 


(2-26) 


Our  main  concern  in  this  section  is  to  derive  a  more  tractable 
radiation  condition,  as  mentioned  before,  which  would  be  a  boundary 
condition  of  Dirichlet,  Neumann  or  mixed  type.  If  one  takes  the 
boundary  to  be  normal  to  the  x-axis,  then  (2-26)  becomes  a  mixed 
type  boundary  condition,  since  in  this  case  we  can  replace  by 


9% 


The  equation  (2-26)  can  he  written  as 


/{«S  ;  _  1)  (  tfm. 

a-A  '  <  f  *  "  < 


(2-26') 


It  is  of  interest  to  note  that,  for  the  infinite-depth 
case,  one  can  obtain  the  coupled  relation  (2-27)  below  for  a 
more  general  boundary  than  (2-26)  or  (2-26').  And  this 
boundary  could  be  any  simply  connected  line  from  the  free  surface 
to  a  sufficient  depth  that  could  represent  the  infinite-depth 
problem  for  numerical  computation: 


(  fT  )  sr  V  C  -7t,  9f  +  TLI  )  , 

(2-27) 

■hiff)  CTt, 
or 

&(<(?>  -»■**$  «-»■*.  r; 

-»-n*9sf  =  vn,f’ 

where  tl  =  (n^ ,n2>  is  the  outward  normal  vector  from  the  fluid 
on  the  boundary. 

It  should  be  noted  that  (2-27)  or  (2-28  )  reduced  to  (2-26) 
or  (2-26')  if  the  boundary  is  chosen  to  be  vertical,  i.e,, 
tl  =  (1,0). 


> 

(2-28  ) 

? 
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2.  Three-dimensional  Problem 

As  mentioned  before,  most  of  th<.  procedures  in  this  section 
are  similar  to  the  previous  section  in  this  chapter.  Thus  a 
brief  treatment  will  be  given.  However,  a  derivation  of  a  new 
radiation  condition  will  be  given  in  some  detail. 

Cylindrical  coordinates  will  be  used  throughout  this  section. 

a)  Local  Disturbance 

The  finite-depth  problem  will  be  treated  first.  It  is 

assumed  that  the  motion  (normal-velocity  distribution)  is  given 
* 

on  R^R0  (refer  to  Fig.  11), 


In  order  to  solve  this  problem  with  a  simple  geometry,  one  can 
use  the  method  of  separation  of  the  variables  as  in  the  previous 
section.  Then  one  obtains  the  elementary  solutions  (Wehausen, 
p.  475)  as  follows: 


*  See  the  footnote  on  page  17 


34 


C^>ik  Tr*(y+H)C/l  +P  YhfaUj  C*S(lU+Tj  t+sfr+Zj' 

cos  -it;  fc+R)[>|  J^HV  t£  KnimO]  <*s(tu+£)  (c*s#n)  <  2"29> 

(i  -  i *2, • • • ) 

where  A,  B,  6,  't  are  constants,  where  the  eigenvalues 
m  ,  m.  (i  =  1,2,’**)  are  defined  in  (2-3)  and  where  n  is  an 

O  1 

integer  because  of  the  juncture  condition  (continuity  condition) 
along  the  (^-direction  unless  there  is  a  radial-plane  cut  pivoted 
on  the  y-axis. 

Without  actually  constructing  a  solution  for  a  given  normal 

velocity  distribution  on  the  body  by  using  the  above  eigenfunctions, 

one  can  examine  the  behavior  of  a  local  disturbance  in  general. 

In  the  first  equation  of  (2-29),  is  the  only  bounded  function 

and  represents  a  standing  wave  with  a  fixed  phase  at  infinity, 

whereas  one  can  find  bounded  standing  waves  of  arbitrary  phase 

in  two  dimensions.  Whenever  one  wants  to  construct  a  wave  of 
* 

arbitrary  phase  at  infinity,  one  must  also  admit  the  singular 
Bessel  function  Y*-  >  in  other  words,  one  should  permit  a  log¬ 
arithmic  singularity,  at  the  axis  ft=0.  Therefore,  the  first 
equation  in  (2-29)  can  be  used  to  construct  a  propagating  wave. 

In  the  second  equation  in  (2-29),  I  is  bounded  at  R=0  but 

’  n 

increases  exponentially  as  whereas  K  is  singular  at  R=0 

n 

but  bounded  for  other  values  of  R  and  decreases  exponentially  as 

*  Only  when  there  is  a  standing  wave  of  arbitrary  phase  can  one 
construct  an  outgoing  wave  which  can  be  expressed  in  the  form 
f -m £<  t-ci:)  for  one  spatial  dimension  and  in  the  form 

<t>  for  two  sPatiLal  dimensions. 
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R  -*»o  .  There  is  no  function  which  is  bounded  everywhere  and 

represents  a  local  disturbance  in  three  dimensions.  In  general, 

one  may  expect  that  the  solution  of  a  boundary -value  problem  will 

require  the  functions  Kn  and  hence  a  singularity  at  R=0.  The 

function  K  will  represent  the  local  disturbance, 
n 

In  order  to  invest  irate  the  behavior  of  a  local  disturbance 
one  may  then  simply  examine  the  function  K^m^R).  The  functions 
( x) ,  and  K^Cx)  have  well  known  asymptotic  representations 


as  x  -*  eo 


.  These  asymptotic  representations  are 


Y-n.  fc)  — '  Srn  (*■-£''*}  j 


(2-30) 


where 


Tilt  TC 

°  T  +  T 


When  n  is  quite  small,  say,  n=0,l,2,3,  these  asympfcotTc 
formulas  give  very  good  approximations  for  these  real  functions, 
even  when  the  variable  X  is  not  so  large.  For  example,  when 
n  =  0  the  asymtotic  formulas  for  and  ^60  give  such  good 

approximate  ones  after  about  one  wave  length. 

When  an  axisymmetric  body  is  heaving  vertically,  n  =  0, 
and  when  the  body  is  surging,  n  =  l,  and  so  on.  The  more  com¬ 
plicated  one  makes  the  motion,  the  higher  the  order  of  the  Bessel 
functions  contained  in  the  solution.  We  restrict  ourselves 
to  the  case  n  =  0  for  simplicity,  but  in  principle  a  general 
type  of  motion  can  be  treated  in  the  same  way  and  there  will  be 
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only  a  change  in  the  region  in  which  the  asymptotic  formulas 

will  be  useful  for  the  approximation. 

From  the  above  considerations  one  can  readily  adopt  the 

two-dimensional  results  (2  —  9)  - (2  “*3)  for  a  local  disturbance 

expressed  in  C  &U  K)  by  replacing  X  by  R. 

For  infinite  depth,  we  shall  examine  a  source  with 

pulsating  strength  cos  Ct  at  y=b  on  the  y-axis 

(Wehausen,  p.  477).  The  potential*  is 

f'i  -V# 

3,l>>  *  tr'4- £  np 

*  (2-31) 

**  tffcci-KO*)  «•*<*]  9 

where 

r1  -  C  +  q-tf  j  nA  *  V  *  (a  V  . 


For  convenience  we  define 

b)  =  +  $  ,  (2-32) 

=<fW  ,  (2_33) 


where 


/-D  f 

ftt  -  iv£yJ  £  . 

For  the  computation  of  the  integral  in  (2-34)  the  Gauss- 
Laguerre  quadrature  formula  has  been  used  as  was  done  for  the 
two-rdimensional  computations  (see  Appendix  A). 


*  If  we  take  a  higher-order  singularity,  then  we  get  higher-order 
Bessel  functions. 


37 


b)  Propagating  Waves  and  Radiation  Condition 

We  have  already  remarked  in  the  previous  sub-section 
that  one  must  have  a  pair  of  and  *n  order  to  construct 
a  p»opagating  wave  and  that  the  asymptotic  formulas  for  , 

Yr,  .  and  Kn  have  a  very  wide  range  of  validity  as  approxi¬ 
mations.  Thus  we  shall  construct  a  new  relation  between  the  , 
approximate  expressions  for  and  ^  ,  which  will  be  used 
as  a  radiation  condition  later.  Since  the  radiation  condition 

is  nothing  but  a  condition  for  the  departing  phase  in  space, 

-1 

one  may  suppress  the  factor  R  in  the  asymptotic  representations 
in  considering  the  oscillatory  parts.  Now  let  us  define  the  first 
term  of  the  asymptotic  expansions  as  functions  to  be  used  in  the 
derivation  of  a  new  form  of  the  radiation  condition  more  suitable 
for  numerical  calculations: 


>  (2-35) 

Y«  ^  o*  -f„>  f 


where 


where 


As  R_>eo 


JnO>R)  -  ZTr\0>RJ  — >0 
Nr,  (  UK)  —  YV,  (ygj  — ^  q 

A  s.  -  s?  . 

Define  tpc  and  as  follows: 

lS  =  Jr*  ■  Jxs  =  A  asCvR-S'n)  y 
T*  "  \J"fc  in  —  A  S"  0^.-9n)  y 


A  ~  p’  = 

"  \Ttp 


constant. 


(2-36) 
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From  the  condition  of  departing  phase  and  (2-36),  in  a  manner 
similar  to  the  derivation  of  the  two-dimensional  case  we  obtain 


C?C)R  -  , 

(<f\  -  Pf  . 


(2-37) 


From  (2-36)  , 

=  5r(J?  a.;=  ±£*3, 

<■9%  .  (2-38> 


From  (2-37)  and  (2-38),  we  obtain  the  following: 


(2-39) 

^  C  -0, 

4T  (  ^)R  +  -V$J  *0.  <2-40> 


We  finally  propose  as  a  new  radiation  condition  for  numerical 
computation 


+  3(t'fp  ~v'fr] 


(2-41) 


where 
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} 


c  =  'fCf  **  ^  (tM.y)  ct , 


At  this  stage  one  may  start  to  ha'. s  doubts  about  the  new  radiation 
condition  (2-41)  because  the  second  terms  (in  the  bracket)  in 
(2-41)  are  0( whereas  only  the  first  term  in  the  asymp¬ 
totic  expansions  for  and  N£t6'R)  was  taken  in  the  derivation 

of  (2-41),  If  the  neglected  terms  are  of  order  less  than  or  equal 
to  0(*"*>  in  the  asymptotic  expansions,  those  terms  should 
have  been  taken  into  account  in  the  derivation  of  (2-41),  In 


order  to  examine  this,  we  write  down  the  complete  asymptotic 
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and  should  have  teen  taken  into  account  in  the  derivation  of 
(2-41)  if  all  terms  through  this  order  are  retained.  Fortunately 
these  worries  turn  out  to  be  unnecessary,  for,  if  one  substitutes 
these  two-term  expansions  (2-44)  into  (2-41),  the  terms  of  0(&  **) 
introduced  by  retaining  the  second  term  of  0(R"^cance*  each 
other,  and  the  lowest  order  left  is  0  .  Thus  we  have 

the  following: 

(x^  +  sjfx*  «  0  cr%)  , 

i  ~  .  (2-45) 

(Y„>«  *  s’rY"  -  >  •' 

As  mentioned  earlier,  the  first  term  in  the  expansions, 

i.e.,  the  equations  (2-35),  are  very  good  approximations  for 

and  »  but  we  take  the  two  terms  as  defined  in  (2-44) 

the  approximations  will  be  far  better.  Even  most  mathematical 

tables  do  not  give  the  Bessel  functions  J  and  Y  for 

o  o 

for  example,  but  give  the  formula  (2-44)  Instead  (e.g.,  McLachlan, 
pp,  215-217),  This  is  even  less  than  three  wave  lengths  in  our 
problem  (x  =  >>R).  Therefore  the  relations  (2-41)  are  valid  to 
the  same  extent  that  (2-44),  neglecting  the  terms  of  Ott 
is  valid.  This  is  then  the  new  radiation  condition  that  we  shall 
use. 

If  we  take  the  limit  (2-41)  reduces  to  the  conventional 
radiation  condition  (1-34-).  Therefore  one  can  interpret  this  new 
radiation  condition  as  a  boundary  condition  which  gives  a  better 
approximation  than  the  conventional  radfatTVn  condition  and  does 
not  require  R  to  be  so  large  for  its  application. 
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If  one  takes  a  boundary  R  =  constant,  the  derivative  with 
respect  to  R  in  (2-41)  becomes  a  Neumann-type  boundary  condition, 
which  we  were  hoping  to  obtain  for  later  application  of  a 
variational  principle. 

It  is  also  of  interest  to  note  that  the  new  radiation 
condition  (2-41)  can  be  obtained  without  using  the  elementary 
•olutions,  i.e.,  the  Bessel  functions  and  ,  but 

through  physical  observations  and  much  simpler  mathematics. 

When  we  observe  steady-state  outgoing  ring  waves  in  a  pond 
which  are  generated  at  a  point  or  in  a  small  region  in  the 
center  of  the  pond,  at  some  distance  away  from  the  wave-making 
region  (or  point)  wo  can  assume 

"T"  **  AdO  j  (2-46) 

where  Y  is  the  elevation  of  the  free  surface  and  the  amplitude 
A  is  a  function  of  R.  Since  we  have  assumed  that  the  motion  of 
the  fluid  is  a  steady  oscillatory  one,  then  the  energy  flux 
across  an  arbitrary  concentric  circular  cylinder  in  one  period 
must  remain  a  constant.  From  this  consideration  we  obtain 

A(R)=?  C/fc  (2-47) 

where  c  is  a  constant. 
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If  we  had  started  with  the  above  equation  (2-47)  In  the 
derivation  of  the  new  radiation  condition,  we  could  have 
obtained  the  sane  equation  (2-41)  without  having  considered 
the  asymtotic  behavior  of  the  Bessel  functions.  This  method 
does  not  seem  to  be  mathematically  rigorous,  but  is  still 
based  on  observation  of  the  real  phenomena,  and  hence  it 
shouldn't  be  disparaged. 
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III.  VARIATIONAL  METHOD 

Most  problems  of  applied  physics  encountered  in  engineering 
are  formulated  in  the  form  of  a  differential  equation  which 
governs  the  behavior  of  a  typical,  infinitesimal  region  in 
the  domain  of  definition,  and  which  can  be  solved  analytically 
or  numerically.  To  solve  such  problems  analytically  one  may 
use  the  method  of  separation  of  variables  for  a  very  simple 
geometry  or  the  method  of  integral  representations  by  using 
Green  functions.  In  s1’ ip-hydrodynamics  problems,  the  latter 
is  used  most  often.  Sometimes,  however,  there  are  difficulties 
in  carrying  through  computations  which  are  Involved  with  Green 
functions.  Furthermore,  the  method  has  another  drawback,  in 
that  it  Is  not  practical  for  a  very  complicated  boundary  geometry, 
i.e.,  for  example,  for  a  variable-depth  problem,  even  though 
it  could  be  done  in  principle. 

In  this  chapter  we  shall  discuss  an  alternative  method, 
a  so-called  variational  method.  A  formulation  in  a  variational 
form  can  be  obtained  directly  from  the  fundamental  physics  of 
the  problem,  e.g.,  the  energy  method  in  a  structural  problem. 

In  a  slightly  different  way,  it  can  also  be  obtained  mathema¬ 
tically  from  the  fundamental  differential  equations.  It  is  not 
always  possible  to  find  a  variational  form  for  a  given  problem. 
When  a  variational  form  is  not  known  for  a  differential  equation 
that  we  wish  to  solve,  we  can  still  use  a  numerical  technique  to 
minimize  a  ’pseudo-variational  form'  or  any  approximate  functional 
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constructed  for  the  differential  equation  (Zienkiewicz,  pp.  38-40, 
1971). 

There  exist  variational  principles  for  the  motion  of  a 
fluid  either  in  an  Eulerian  representation  or  in  a  Lagrangean 
representation  (Eckart,  1960;  Luke,  1966).  Since  our  formulation. 
Equations  (1-18)  -  (1-23)  and  the  radiation  conditions  (2-26*) 
and  (2-41),  consists  of  Laplace's  equation  with  boundary  conditions, 
partly,  of  Neumann  type  and,  partly  of  mixed  type,  we  shall 
simply  adopt  a  well-known  variational  form  for  this  problem 
(Mikhlin,  pp.  138-151,  1964). 

Let  us  consider  Laplace's  equation 

V*  $  CX.*))  —  0  Ct*  (3-1) 

with  a  boundary  condition,  of  the  form 

<j)^  •+  oi  ^  «  p  C*T»  >  j 

where  o((x>^)  and  are  known  functions.  The  variational 

form  for  this  problem  requires  introduction  of  the  following 
functional: 

FM>)  =jf  +|c*w<P  -f  W* .  <3-2) 

A  function  that  minimizes  the  functional  F  is  a  solution 

of  (3-1),  and  vice  versa.  The  proof  that  the  function  which 
minimizes  the  functional  is  the  solution  of  (3-1),  is 

also  in  Mikhlin  (pp.  115-116,  1964). 

In  three  dimensions  one  has  a  similar  form.  Let  satisfy 


(in  SL  ) , 
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(on  aft. )  f 


(3-3) 


where  oifiyjjjl  and  are  known  functions.  Then  the  corresponding 

function*)  H4>)  is 


F(4>)  =J[  i-KridjJ, 


(3-4) 


It  is  of  interest  to  note  a  special  case  in  which  the 
general  form  of  the  functional  in  (3-2)  or  (3-4)  can 
degenerate  to  the  boundary  integral  only  when  we  distribute 
Green  functions  along  the  boundary  (Bessho,  1970),  The 
application  of  this  degenerate  case  is  given  in  Sao,  Maeda  and 
Hweng  (1971).  In  this  case  the  computation  will  be  similar  to 
the  scheme  used  by  Frank  (1967), 

In  our  problem,  from  (1-18)  -  (1-23)  we  can  readily  see  that 
we  will  have  two  functionals,  one  for  a  cosine-mode  and  the 
other  for  a  sine-mode,  and  that  these  two  functionals  to  be 
minimized  are  coupled  on  the  boundary  on  which  the  radiation 
condition,  (2-26')  or  (2-41),  is  imposed,  i.e.,  the  integral 
along  the  boundary  in  (3-2).  Due  to  this  coupling  on  a  boundary, 
the  final  matrix,  which  is  the  coefficient  matrix  in  a  set  of 
algebraic  equations  obtained  from  the  finite-element  discreti¬ 
zation  to  be  discussed  in  the  next  chapter,  is  not  symmetric, 
when  we  combine  two  sets  of  the  algebraic  equations,  for  the 
cosine-mrJe  and  the  sine-mode,  respectively,  whereas  in  most 


structural  problems  they  are  symmetric.  However,  we  can  obtain 
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the  two  functionals  f  \  fe)  ,  p S(<f *)  ,  for  the  functions 

?C  and  f'S  ,  respectively,  from  (3-2),  (1-8)  - 
(1-23),  and  (2-26')  and  Fig.  12  as  follows: 


f%c) 


(3-5) 


si 


-f(W V*)  +)>\  ?Vii 


(3-6) 


where  the  Integrals  along  the  bottom  $  vanish  in  both  cases 
unless  the  bottom  is  moving.  As  mentioned  before,  the  last 
integrals  along  the  ' radiation  boundary'  in  (3-5)  and  (3-6) 
show  the  coupling  relation  between  pc  and  therefore 

the  minimizing  functions  and  must  be  found  at  the  same 

time,  not  one  after  another,  unless  we  can  decouple  them. 

C  8 

When  1  pressure  distribution  p  cos  6"  t  +  p  sin  (T  t 
is  specified  in  a  segment  on  the  free  surface,  then  we  obtain 
the  new  functionals  after  we  add 

j  it 

dif 


(3-7) 


to  (3-5)  and  add 


4  (j-V/e (3'8) 

b, 

to  (3-6).  When  we  have  stratified  fluids,  the  boundary  Integrals 

along  the  interface,  for  the  upper  fluid  and  for  lower  fluid 

(see  (1-28)  and  (i-29),  will  show  the  coupling  relation  between 

the  upper  fluid  and  the  lower  fluid  in  each  cosine-mode  and  sine- 
c  s 

mode  functionals,  F  and  F  ,  respectively. 


In  three  dimensions  heaving  motion  of  an  axi-symraetric 
body,  which  can  be  formally  reduced  to  the  two-dimensional 
problem,  will  be  treated  here  for  simplicity.  But  a  general 
three-dimensional  problem  can  be  described  in  a  straightforward 
way.  From  the  formulations  in  the  first  chapter  and  the  new 
radiation  condition  (2-41),  we  obtain  two  coupled  functionals 
and  F *(.*$*)  (refer  to  Fig.  1.3)  as  follows: 


'c?‘>  =  t-  - 


-2*^  -21 , 
w  v>-  «$ ^  ^ P.Crt'Jf 

-Ji (t,(n -v!) +«j.(yVfs<ij  >  <3 

or,  by  dividing  (3-9),  (3-10)  by  1C  and  redefining  F°  and  FS, 


(3-9) 


fc«') = #  tG^+^jyR  -  ttffi*.  <3-u) 

FS(<fs)  -  R  [l$‘  -K^J  V*  - 


(3-12) 
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IV.  NUMERICAL  PROCEDURES 

In  this  chapter  a  brief  description  of  the  method  of 
finite-element  discretization  will  be  given  in  the  first 
section.  How  an  infinite  boundary  can  be  truncated  in  the 
numerical  procedure  will  be  discussed  in  the  second  section. 

1.  Finite-element  Discretization 

A  reader  can  find  a  very  extensive  and  detailed  exposition 
in  Zionkiewicz  (1971),  Therefore  we  will  give  only  a  very 
brief  description  here.  The  convergence  of  the  numerical 
solution  of  Laplace's  equation  is  discussed  by  Reid  (1972). 

To  begin  with,  the  functional  F(<^)  of  a  single  function 
defined  in  (3-2)  will  be  treated  first  with  the  goal 
of  developing  a  numerical  procedure  for  finding  the  function 
<|>  that  minimizes  the  functional  F  .  The  coupled  case, 
which  we  have  in  (3-5)  and  (3-6)  or  (3-11)  and  (3-12),  will 
be  discussed  later. 

Let  the  region  occupied  by  fluid,  up  to  the  place  at  which 
the  radiation  condition  is  to  be  imposed,  be  subdivided  by  lines 
or  surfaces  into  a  (not  necessarily  rectangular)  grid.  Each 
connected  piece  within  the  subdivision  will  be  called  an  'element'. 
We  suppose  to  be  a  function  that  is  continuous  and  bounded 
(but  see  the  discussion  in  chapter  VII)  in  the  subdivided  region. 

*  The  function  4*  defined  in  (3-1)  and  (3-2)  can  be  any 
integrable  function  and  it  should  not  be  understood  as  the 
velocity  potential  defined  in  (1-33), 
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One  of  the  important  steps  in  the  procedure  is  the  introduction 
of  a  set  of  interpolation  functions  N^(x,y),  i  =>  1,»*», 

associated  with  each  element  and  of  such  a  character  that  can 
be  approximated  as  a  sum  of  these  functions,  each  multiplied 
by  the  value  of  at,  say,  a  node  of  the  grid  associated  with 
the  element  (  ;  at  the  i-th  node).  However,  these  values  of 

need  not  be  nodal  values  of  but  nay  be  other 

values  (parameters)  characterizing  in  the  element,  for 
our  numerical  scheme  requires  minimizing  a  functional  which 
is  represented  in  integral  form  rather  than  directly  with  the 
values  of  function  itself.  Let  us  write  the  set  of  interpolation 
functions  as  a  row  vector 


[u]  =[Ni,Nx,  ♦  •  •  ,  Mm] 


(4-1) 


and  the  set  of  'nodal'  values  as  a  column  vector 


itf 


•  ( $T,  4,  •  •  (4-2> 

in  an  N-node  element.  The  superscript  e  on  W  or  on 
<^(  means  that  these  values  are  considered  in  an 

individual  element.  We  may  then  approximate  <|>  in  each 
element  by  the  sum 


(4-3) 


We  shall  give  below  an  example  for  a  rectangular  element, 
a  very  simple  element  shape,  for  the  purpose  of  illustration. 
However,  in  our  actual  computations  we  used  more  elaborate  ele¬ 
ments,  a  four-node  quadrilateral  and  eight-node  quadrilateral 
(see  Appendix  B). 
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1-6,' 

J 

(«,k) 

i 

4 

2. 

_ 1 

( 

l 

We  can  write  an  In  (4-3) 

<j>  5*  Mi  $  +^4..  +  ^  4i  ^  ^  4  }  (4-4) 

where 

M,  -=  «•- >  , 

kL  (XV  ~  ~~  («.-»  (t-j)  , 

.  (4-5) 

Nj  fry  ®  qtfc  <«**)  (H)  , 

(*|J)  =  (CbfX)  fap  # 


In  order  to  minimize  the  functional  F(<^)  1**  (3-2)  with 
respect  to  the  total  number  of  parameters  (or  nodal  values) 
associated  with  the  whole  domain,  we  can  write  a  system  of 
equations 


9F 

»|+f 


V 


/ 


0 


(4-6) 
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Let  an  element  in  Sit.  be  denoted  by  .  Then  we 

decompose  JZ  as  follows; 

F»'-l  *L  -fit  *JJ 

=  1 

e 

We  now  approximate  within  each  element: 

f*<+)c  e*c 


(4-7) 


(4-8) 


Henceforth  we  shall  simply  write  F  for  the  approximate  value. 
It  now  becomes  evident  that  the  interpolation  functions 
must  be  chosen  in  conformity  with  the  nature  of  the  functional 
F  ,  In  particular,  they  should  be  chosen  so  that  at  the 
interfaces  the  approximation  to  ^  is  such  that  it  and  Its 
derivatives  of  one  order  less  than  those  occurring  In  the  inte¬ 
grands  of  p  are  continuous.  This  assures  that  there  Is  no 
contribution  to  the  Integral  from  the  Interfaces. 

Then  from  (4-6)  and  (4-7)  we  obtain  (4-9) 


-2E-«  iff  = 0 


(4-9) 


For  any  node  we  can  write,  by  differentiating  (3-2)  with  respect 
to  ^  (1=1,2,...), 

ifr  =|  *♦•&*■*» 

+j  (4-l0) 
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where  the  second  integral  is  present  only  if  the  element  has 
a  boundary  on  which  the  boundary  condition  in  (3-1)  is 
specified. 


Noting  that  is  no  longer  a  function  of  x  and  y 

but  that  [N]  is  now  a  function  of  x  ands  y, 


it 


f 


(4-11) 


finally  we  obtain  in  the  whole  region 


-M=o. 

041  $  =[b1  , 


where 


f 


(4-12) 


(4-12’ ) 


(4-13) 


(4-14) 


Mt  «ls  . 


(4-15) 


The  formulated  problem  has  now  been  reduced  to  the  equation 
(4-12')  a  set  of  linear  simultaneous  algebraic  equations. 

The  coefficient  matrix  [A]  has  the  nice  properties  of  being 
symmetric,  as  one  can  see  readily  from  the  equations  (4-13) 
and  (4-14)  and  of  being  banded  if  nodes  are  properly  numbered. 

The  numerical  computations  of  the  integrals  (4-13),  (4-14) 
and  (4-15)  for  four-  and  eight-node  quadrilateral  elements  are 
briefly  discussed  in  Appendix  B, 

Now  let  us  consider  the  finite-element  discretization 
for  the  coupled  functionals  JH *  y  ^  and  in  (3-5) 

and  (3-6)  or  (3-11)  and  (3-12),  Since  we  can  use  (4-8)  In 
two  dimensions  without  change  and  also  since  we  can  write  with 
a  slight  change  in  three  dimensions 


hi 


(4-16) 


where  N^(R,y),  let  us  consider  next  the  boundary  integral 

involving  ,  defined  in  (3-1)  and  reduced  to  (4-14)  in 

the  final  matrix  equation  (4-12').  In  two  dimensions  ^  occurs 
only  on  the  free  surface,  where  o(=s-®4t  ,  a  constant  .  Then, 


*  When  we  have  two  fluids,  we  also  have  at  the  interface  of  the 
two  fluids,  ©(.  and  (see  (1-28)  and  (1-29), 
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referring  to  Fig.  12,  we  obtain  from  (4-14) 


hi 


ix 


(4-17) 


both  for  *=c(<f‘)  and  F  (?’>  . 

In  three  dimensions  there  are  two  boundaries  on  which 
occurs,  on  the  free  surface  and  on  the  'radiation  boundary'. 

Q 

Thus  we  obtain  an  analogous  form  to  (4-14)  for  both  F  and 
FS  (see  Fig.  13): 

S'  *  4  * • 


(4-18) 


Next  let  us  consider  the  boundary  integrals  involving 
P  (x/y).  The  value  of  <4  considered  above  was  constant  in  our 
problem,  but  |3  may  be  a  function.  We  can  express  analogously 


to  (4-3)  as 


p-W  {f} 


V 


where  [N]  are  the  interpolation  functions  as  before  and 
are  the  nodal  values  of  the  function  jl  .  From  (4-15)  and 
(4-19)  we  obtain 


(4-19) 


k  =  I  (Ml  [f}  MiJs  =  (j&O  •  <4-20) 


The  equation  (4-20)  becomes  in  two  dimensions 


b;  =  ( 


CtJI  Niil)’  fj  . 


!  i 


i  I 
! 


3 

* 

% 


l 

£ 

% 


(4-21) 
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Vi. 


The  values  of  jflj  ,  that  is,  the  normal-velocity  distributions 
are  known  and  the  computation  of  the  integrals  along  the 
body  boundary  in  (4-20)  and  (4-21)  is  straightforward.  These 
become  the  components  of  [B]  in  (4-12').  For  convenience  let 
the  matrix  equation  be  understood  such  that  the  right-hand 
side,  [B],  is  known,  whereas  the  unknown  values  are  on  the 
left-hand  side,  as  is  conventional.  Now  let  us  consider  the 
integrals  along  the  boundary  ^1  +^2  or^  in  the  expressions 
(4-22)  to  (4-25)  for  both  cases,  i.e.,  the  cosine  mode  and  sine 
mode,  the  boundary  integral  contains  the  values  of  the  sine-mode 
potentials.  First  we  compute  all  and  next  we  move  those 
terms  which  contain  the  opposite-mode  potentials  to  the  left-side 
hand.  In  order  to  combine  the  cosine-mode  potentials  and 
the  sine-mode  potentials  JP*  we  rearrange  both  potentials  into 
one  array  as  follows: 
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(4-26) 


This  arrangement  minimizes  the  bandwidth  of  the  coefficient 
matrix  after  combination.  The  coefficient  matrix  will  now  be  no 
longer  symmetric. 

Further,  if  we  denote  the  integrals  along  and  In  (4-22) 
and  (4-25)  as  b^0  ,  and  those  in  (4-23)  and  (4-25)  as  b^S,  the 
final  form  of  the  matrix  equation  becomes: 


(4-27) 


When  an  oscillatory-pressure  distribution  is  specified 
on  a  segment  of  the  free  surface  or  when  there  are  two  fluids 
of  densities  f,  ,  ,  all  the  procedures  for  obtaining  (4-27) 

will  be  very  similar  to  those  explained  earlier.  It  is  of  interest 
to  note  that  the  interface  condition  in  the  two  fluids  (see  (1-28) 
and  (1-29))  makes  the  coefficient  matrix  [A]  in  (4-27)  asymmetric 
dije  to  the  coupled  interface  boundary  condition  between  the 
upper  fluid  and  the  lower  fluid. 


2,  Truncation  of  the  Infinite  Boundaries 


In  Chapter  II,  the  decaying  behavior  of  the  local  disturbance 
has  been  considered  and  a  new  radiation  condit Lon> was  derived 
which  can  be  applied  where  the  local  disturbance  is  negligibly 
small  compared  with  the  propagating  waves.  We  can  state  a  few 
criteria  for  the  subdivision  of  the  domain  occupied  by  the  fluid 
(or  fluids)  and  for  the  truncation  of  the  infinite  boundaries. 
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Criterion  1,  One  should  subdivide  each  wave  length  by  at 

least  ten  approximately  equidistant  points  along 
the  horizontal  direction. 

Criterion  2,  The  truncation  of  the  infinite  boundaries  should 
be  made  after  examination  of  the  behavior  of  the 
local  disturbance  discussed  in  Chapter  II. 

It  is  difficult  to  state  definite  and  detailed  criteria 
for  Ihe  general  problems,  but  one  can  sensibly  determine  the 
truncation  of  the  Infinite  boundaries  and  the  proper  subdivisions 
in  the  domain  case  by  case  after  understanding  the  specific 
problem  which  is  dealt  with.  For  example,  if  the  bottom  is  not 
uniform  as  shown  in  Fig.  14,  then  the  subdivisions  and  the 
truncation  of  the  infinite  boundaries  should  be  different  along 
the  left-hand  side  and  the  right-hand  side  as  shown  in  Fig.  14, 
since  we  know  the  asymptotic  form  of  the  propagating  waves  on 


FIG.  14.  Subdividing  Meshes  in  the  Fluid 


In  this  case  )?  in  the  radiation  condition  is,  of  course, 
different  on  the  left-hand  side  and  on  the  right-hand  side, 
i.e.  ynfy,  on  the  left-hand  side  and  \>-*Vx  on  the  right- 


hand  side,  where 
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V.  TESTING  OF  THE  METHOD 

1.  Pure  Propagating  Waves 

In  this  section  simple  tests  are  made  to  examine  the 
radiation  condition.  Let  us  consider  a  fictitious  wall  on 
which  a  normal  velocity  is  defined  as 

mv  «  ($P)„  (5-D 

where  <|>p  is  defined  in  (2-19)  in  two  dimensions. 

The  numerical  results  for  the  boundary  condition  (5-1)  on 
;t  fictitious  wall  have  been  compared  with  the  exact  solution 
(2-19)  for  many  cases.  For  example,  in  one  case  the  radiation 
condition  is  applied  right  next  to  the  wall;  in  other  words, 
one  discretized  element  has  both  boundaries,  a  fictitious- 
motion  boundary  on  one  side  and  the  radiation  boundary 
condition  on  the  opposite  side.  All  the  numerical  results 
give  very  good  agreement  with  the  exact  solutions  whenever  we 
subdivide  the  domain  of  the  fluid  properly  into  finite  ele¬ 
ments.  Consider  a  thin  strip  one  hundredth  of  a  wave  length 
wide  subdivided  into  one  hundred  elements  in  the  vertical 
direction.  If  the  fictitious  motion  is  applied  at  the  left- 
hand  side  of  the  strip  and  the  radiation  condition  at  the 
right-hand  side,  then  the  numerical  results  agree  with  the 
exact  values  to  six  or  seven  decimal  places  if  eight-node 
quadrilateral  elements  are  used. 

In  three  dimensions,  propagating  ring  waves  are  tested 
with  the  local  disturbance  artificially  suppressed, 
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i.e.  the  boundary  condition  on  a  fictitious  vertical 
cylinder  R  =  R  is  taken  as  follows: 


“V 


c^n  vH 


J.  (vO  , 


3fL|e.R<>- 


YCvO  . 

c^sh  yH 


(5-2) 


where  H  is  the  depth. 


The  exact  solution  for  the  boundary  condition  (  5-2) 
is  simply 


?  <■  R.'j.'t)  =*  Cl?  +  «■*»• 


— — -  V  ^ —  (x  Os6t  -vX  (5-3) 

co* h  UH 


where  R  2R#  -H  £  ^  i.0  ,  In  this  three-dimensional 

case  the  numerical  results  are  compared  with  the  exact  solu¬ 
tions  (Bessel's  functions)  given  by  Abramowitz  and  Stegun  (1967) 
and  shown  in  Fig.  20  for  y=0.  The  radiation  condition  is 
applied  at  RsR^,  which  was  taken  as  less  than  one  and  a  half 
wave  lengths  for  the  two  cases  considered.  The  numerical 
results  are  identical  up  to  three  significant  digits  with  the 
Bessel  functions. 

2.  Pure  Local  Disturbance 

As  we  have  seen  in  Chapter  two,  if  we  apply  the  radiation 
condition  on  a  boundary  which  is  so  close  to  the  moving 
body  that  the  local-disturbance  potential  <5.  has  not 


decayed  sufficiently,  then  we  will  obtain  a  numerical 
solution  not  for  the  original  problem  but  for  a  completely 


different  problem.  For,  when  we  take  the  conventional 
radiation  condition,  which  states  the  asymptotic  behavior 
at  infinity,  then  the  local  disturbance  potential 

trivially  satisfies  the  radiation  condition. 

Therefore  no  special  difficulties  in  connection  with  the  local 
disturbance  arise  in  an  analytic  method.  On  the  other  hand, 
in  numerical  scheme,  the  treatment  of  the  local  disturbance 
is  the  most  cumbersome. 

In  this  section  the  radiation  condition  is  applied 
at  various  distances  at  which  the  local  '•isturbance  remains 
of  considerable  magnitude,  in  order  to  obtain  some  idea  as  to 
how  close  the  truncation  of  the  infinite  boundary  may  be  applied. 

Let  us  consider  an  example  in  two  dimensions  with  the 
exact  solution 

H  €  *'*  Hf»i  C9*H)  LosSt  ^  (5-4) 


where  PI,  is  given  in  (2-3)  and  A  is  a  constant, 
'  o 

The  two  components  of  £  are 


?c<^y)  =  K  *  m,x  <**  , 

f 1 =  o  , 


(5-5) 


where  -HiyiO  and  04x<*°.  Then  on  a  fictitious  wall  at 
x=0  the  boundary  condition  will  be: 
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fy  ly,0  “  ~m'  ^  W.  , 

£l*«  =  «■ 


(5-6) 


We  apply  the  radiation  condition  at  x=x^ . 

Fig.  21  shows  the  potentials  *f  and  with  the 
radiation  condition  imposed  at  various  distances  from  the 
plane  x-0  and  shows  that  the  numerical  solution  agrees 
very  well  in  this  case  when  the  radiation  boundary  is 
taken  such  that  x/H£3.  However  when  we  impose  the  radiation 
condition  too  close  to  the  moving  body,  then  the  numerical 
solution  gives  a  non-zero  sine-mode  potential 
although  this  is  identically  zero  in  the  exact  solution. 


3.  Oscillatory  Pressure  Given  on  the  Free  Surface 

Let  an  oscillatory  pressure  be  given  on  a  segment  zt 

of  the  free  surface;  the  boundary  condition  on  the  free 

surface  is  then  given  in  (1-31).  Stoker  (PP.  58-66)  has 

applied  complex-variable  theory  to  give  the  solution  for  the 

propagating  waves,  i.e.  for  25  ,  as  an  asymptotic  solution 

r 

at  a  large  distance.  If  we  adopt  the  Green  function  for  this 
problem  and  make  use  of  Green's  theorem,  the  solution  can 


be  given  as 


(« 


(5-7) 
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where  G  is  defined  in  (2-14)  and  C  is  defined  in  (1-31). 

For  finite  depth  G  is  given  in  Wehausen  (1960);  it  would 
not  be  difficult  to  extend  (5-7)  to  three  dimensions. 

In  this  section,  two  cases  in  water  of  infinite  depth 
are  compared  with  the  solution  computed  from  (5-7)  and  given 
in  Fig.  22.  Fig.  23  shows  the  same  two  cases  for  finite 
depth. 

It  is  of  interest  to  remark  that  a  problem  with  the 
non-homogen#ous  free-surface  boundary  condition  must  be 
solved  whenever  we  treat  higher-order  problems. 

4.  Forced  Motion  of  Two-dimensional  Cylinders  in  Water  of 
Infinite  and  Finite  Depths 

A  circular  cylinder  of  radius  a  oscillating  in  a  free 
surface  has  been  tested  for  infinite  depth  for  five  non-dimensional 
wave  numbers  ya  =  0.5,  1.0,  1.3,  2.0,  3.0,  and  compared  with 
Porter's  results  (Porter,  1960).  They  give  very  good  agreement 
in  all  cases.  Table  1  shows  the  comparisons  between  Porter's 
results  and  the  results  by  this  numerical  method  for  the  case 
ya  =  l. 

For  finite  depth,  1!  =  2a  a  hal f-immersed  circular  cylinder 

has  been  tested  for  two  non-dimensional  wave  numbers 
V  0s-  0.1,  0,5  and  compared  with  the  results  obtained  by  another 
method  (C.  H.  Kim,  1969).  The  added-mass  coefficients  and  wave- 
ampl.' tude  ratio  for  heave  do  not  give  good  agreement  for  either 


case. 


Damping  Coeff.  Asympt.  Wave  Amp.  (Y) ,  Force  (F) 


Table 
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5.  Two-dimensional  Diffraction  Problems  In  Water  of  Infinite 
and  Finite  Depth  (Circular,  Rectangular  Section*) 


The  incoming  wave  is  assumed  tc  be 


Yi 

as  yx  c.*rt:-osvx 

with  the  amplitude  being  unity.  Following  the  definition 
(1-39),  for  the  incoming  wave,  we  write 

Yt  •  V*  «■* c-t  't*  *&**fe* 


(5-8) 


(5-9) 


where 


Yx  ®  *»*»*>*  ,  Ys 

From  (1-40)  we  define  the  free-surface  wave  profile  due  to 
a  fictitious  forced  .notion  given  in  (1-25): 

Yd  cx.tf  *  YoC  «•*  «*+ V  > 


where 


Y>c  =  -  -  j-  Vj  (x,<>)  # 
Yo5-T^  ^  > 


(5-10) 


then  the  total  velocity  potential  ^  is  given  as 

i§r  ~  3k  '*'%> 

and  the  wave  profile  on  the  free  surface  is 
Y  C  X  *)  -  Y  (-X&  -t  Y  Ot.-fc) 

and  v  i 

Yrc*.-t)  =  yrc  t-sr+  +  yT 

In  order  to  describe  the  phase  relation  between  Y,C  and 


(5-11) 


(5-12) 


> 


6B 

we  define  the  phase  angle 

&pfc)—  ' ba-vT1  C  ^ / Y p)  j  (5 

similarly 

&T&)~  W1  C-Vt/Yt)  (g 

where  the  arc  tangent  is  understood  to  take  ita  principal 
value,  i.e.,  S  dp,  Op  -  % 

«  M.tn.  M  are  the  amplitudes  at  infinity  of  the 
Incoming,  transmitted,  and  reflected  waves,  respectively, 
then  we  define  the  reflection  and  transmission  coefficients  by 

cR=iYi/|vT|  _  4=|yti/iy,i  (5 


and  obtain  as  a  consequence  of  conservation  of  energy, 


*■ 

M" 


(5 


we  can  also  express  the  transmitted  wave  at  infinity  similarly 

'"<5-8>.  yr  =  |y|  Sr« 


where  is  the  phase  lag  of  the  transmitted  wave  to  the 

incoming  wave  at  infinity. 

Diffraction  of  homogeneous  incoming  waves  on  a  submerged 
ular  cylinder  was  taken  as  a  test  computation  since  this 
problem  has  been  treated  analytically  by  various  people 
(Dean,  1948;  U’sell,  1950;  Ogilvie,  1963),  It  is  well-known 
that  the  reflection  coefficient  of  a  submerged  circular 


13) 


14) 


15) 


16) 


circ- 


cylinder  in  water  of  infinite  depth  is  zero  in  the  first-order 
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theory.  Thus  the  model  is  taken  for  the  numerical  calculation 
to  be  identical  with  the  one  Dean  and  Ursell  computed,  as  shown 
in  Fig.  15. 

The  results  fory*.=^are  shown  in  Fig.  24.  The 
transmitted  wave  lags  behind  the  incoming  wave  by  about  82 
degrees;  Ursell  found  about  90  degrees. 

A  diffraction  problem  with  a  rectangular  cylinder  (refer  to 
Fig.  16)  in  water  of  finite  depth  has  been  tested  for  some 
frequencies  and  compared  with  the  results  obtained  by  another 
method  (Mei  and  Black,  1969),  Hie  results  show  very  good  agree¬ 
ments  for  both  submerged  and  surface  obstacles.  Fig.  25  -  Fig.  28 
show  the  results  for  these  obstructions. 

Each  figure  giving  the  results  of  the  diffraction  problems 
consists  mostly  of  four  sheets  of  figures;  for  example,  Fig.  25 
consists  of  Fig.  25a,  Fig,  25b,  Fig,  25c,  and  Fig.  25d.  The 
first  figure  shows  the  diffracted  waves  ^W.  The  upper  part 
of  tiie  second  one  shows  the  phase  angle  between  the 

cosine-mode  waveband  the  sine-mode  wave  'jJ1  .  This  figure 
shows  the  transmitted  waves  on  the  right-hand  side  and  the  sum 
of  the  incoming  wave  and  the  reflected  wave  on  the  left-hand  side. 
The  fourth  figure  is  arranged  in  a  similar  fashion  to  the  second 
figure  except  for  the  total  waves.  The  phase  angle  0T  shows 
the  phase  lag  of  the  transmitted  waves  with  respect  to  the 
incoming  waves  at  the  far  right-hand  side;  it  is  hard  to  give  a 
physical  interpretation  to  $j-  on  the  left-hand  side  when  the 
reflection  coefficient  £*,  is  large,  as,  e.g. ,  in  Fig.  32d.  The 
lower  part  of  the  fourth  figure  shows  the  reflection  and 
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transmission  coefficients  at  the  far  left-hand  and  far  right-hand 
sides,  respectively.  It  is  easy  to  prove  analytically  that 
w  max  on  the  left-hand  side  has  a  sharp  trough  and  a  smooth 
crest  and  that  the  average  of  two  values  of  ax  at  the  crest 
and  trough  is  one.  But  sometimes,  for  example,  Fig.  26d  does  not 
show  a  very  sharp  trough,  but  is  chopped  off  in  the  plotting 
by  the  computer  due  to  the  lack  of  a  sufficient  number  of  data 
points  near  the  trough. 
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A  hemisphere  in  heave  motion  in  water  of  infinite 
depth  has  been  tested  for  two  frequencies,  0.5,  1.0, 

and  compared  with  the  results  obtained  by  other  methods 
(Cumming,  1963),  The  results  show  very  good  agreement. 

The  waves  and  wave  forces  generated  by  vertical  circular 
cylinder  of  draft/radius  =  0,5  in  heave  motion  in  water  of 
infinite  depth  have  been  computed  for  few  frequencies,)/^  0.5, 
1.0  ,  The  resu? ts  have  been  compared  with  those  of  Sao,  Maeda 
and  riwang  (l&'/l),  but  they  are  not  in  good  agreement.  In  view 
of  the  good  agreement  between  our  results  and  those  computed, 
by  many  others  for  other  configurations,  we  are  inclined  to 
favo,  our  results  in  case  of  disagreement. 


VI .  NEW  PROBLEMS 


Once  this  method  has  been  proved  to  be  a  useful  numerical 
method,  one  can  utilize  it  for  solving  any  problem  within  the 
scope  of  its  applicability,  however  complicated  its  boundary 
may  be.  In  this  chapter  only  a  few  sample  problems  will  be 
treated  since  any  one  can  make  use  of  the  method  for  cases  of 
interest  to  him. 

1,  Forced  Motion  of  Two-dimensional  Cylinders  in  Water  of 
Variable  Depth 

Forced  motion  of  a  rectangular  cylinder  in  the  free 
surface  with  a  vertical  cliff  submerged  under  the  heaving 
cylinder  is  treated.  As  shown  in  Fig.  17,  we  have  an 
infinite  depth  on  the  left-hand  side  and  a  finite  depth  on 
the  right-hand  side.  Hence  we  have  two  different  lengths  of 
the  waves  which  propagate  to  each  direction.  Due  to  the 
asymmetry  of  the  bottom  with  respect  to  the  y-axis,  we  may 
expect  a  non-zero  force  component  in  the  x-direction,  even 
though  the  cylinder  and  the  motion  are  symmetric  with  respect 
to  the  y-axis.  The  numerical  computation  for 
shows  that  the  amplitude  of  the  x-component  of  the  hydrodynamic 
force  is  about  one  fifth  of  that  of  the  y-component.  Fig.  29 
shows  the  wave  profiles  for  this  case. 
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2,  Diffraction  Problems  In  Water  of  Finite  Depth  (Triangular, 
and  Sinusoidal  Section  and  Vertical  Multi-barrier  along 
the  Free  Surface) 

The  geometries  of  the  triangular-  and  sinusoidal-shaped 
obstructions  on  the  bottom  are  shown  in  Fig.  18.  The  results 
for  the  triangular-shaped  obstacle  on  the  bottom  are  given 
for  yh  =  1.0  in  Fig,  30.  The  results  for  the  sinusoidal  hump 
on  the  bottom  are  given  in  Fig.  31  for  yfi=  0.5. 

A  diffraction  problem  for  two  vertical  flat  platds 
piercing  the  free  surface  was  treated  in  some  detail  and 
this  problem  was  extended  to  the  case  of  two  fluids  of 
different  densities,  and  (corresponding  to  oil 
contained  between  two  vertical  flat  plates  (see  Fig.  19). 
Experimental  results  (Flaissi,  1972)  are  available  for  these 
problems.  The  results  are  obtained  for  four  different 
frequencies  for  a  homogeneous  fluid  and  are  compared  with 
the  experiments.  Agreement  is  good.  The  results  for  two 
fluids  are  obtained  for  a  frequency  which  is  near  one  of 
the  resonant  frequencies  of  the  internal  waves.  The  maximum 
amplitude  of  the  internal  waves  is  about  six  times  higher 
than  that  of  the  surface  waves  on  the  oil.  Fig.  32  shows 
the  results  for  a  homogeneous  fluid  for  period  T  =  1  sec. 

(  yh  *r  0,6404)  and  a/h  a  2.34,  Fig.  33  -  Fig,  35  show  the 
results  for  the  two  fluids  for  T  =  1.84  sec,  (  h=0,2704)  and 


a/h  =  2.34. 


VII.  DISCUSSION 


The  main  purpose  of  this  paper  is  rather  the  testing 
of  a  numerical  method  than  anything  else.  Therefore  the 
testing  of  this  method  has  been  carried  out  for  a  variety 
of  different  problems,  varying  from  homogeneous  wave’ 
propagation  to  a  diffraction  problem  due  to  two  vertical 
flat  plates  piercing  the  free  surface  when  oil  is  contained 
between  two  vertical  plates. 

Most  results  obtained  by  this  method  agree  with  those 
obtained  by  other  methods,  either  analytic,  numerical  or 
experimental.  But  there  were  two  types  of  problems  which 
did  not  agree  very  well:  one  is  the  forced  motion  of  a 
two-dimensional  circular  cylinder  in  water  of  finite  depth 
and  the  other  is  the  heaving  vertical  circular  cylinder  in 
water  of  infinite  depth.  It  seems  to  be  worthwhile  to  compute 
the  complete  results  since  we  checked  only  a  few  frequencies. 
One  can  also  easily  compute  the  forced  motion  of  the  axi- 
symmetric  vertical  cylinder  in  water  of  finite  depth. 

In  addition  to  the  fact  that  any  degree  of  complicated 
boundary  geometry  can  be  handled  easily  by  this  method,  there 
is  still  a  further  degree  of  freedom  in  that  one  can  express 
the  velocity  potentials  on  a  mathematically  singular  plate 
(in  which  the  potentials  are  discontinuous  across  the  boundary) 
through  numbering  of  the  nodes  for  each  element.  One  can 
easily  construct  a  'numerical  Riemann-«urface  branch  cut'  by 
making  two  adjacent  elements  which  have  the  singular  boundary 


n 


in  common  have  different  nodal  numbers  at  that  boundary,  even 
though  the  coordinates  of  the  nodes  on  the  boundary  are  identical. 

At  this  stage  we  may  ask  ourselves,  what  happens  at  the 
ends  of  the  flat  plate  which  are  submerged.  Does  this 
method  give  the  correct  behavior  at  the  singularity?  The 
solution  of  the  velocity  potential  at  the  singular  point 
breaks  down  unless  we  introduce  a  proper  interpolation  function 
to  represent  the  near  field  including  the  singularity.  If  we 
use  a  linear-  or  quadratic-interpolation  function  in  the  near 
field,  then  we  always  obtain  a  finite  value  at  the  singularity, 
because  the  singularity,  which  is  Q(t^) in  the  near  field,  is 
integrable  and  our  numerical  method  minimizes  a  functional 
whi  is  represented  in  integral  form  rather  than  dealing  with 
the  function  itself. 

It  would  seem  to  be  worthwhile  to  develop  a  proper  inter¬ 
polation  function  to  represent  the  near  field,  Including  the 
singularity,  in  order  to  obtain  a  correc*  solution  at  the 
singularity. 

It  would  be  very  interesting  to  extend  this  method  to 
the  ship-resistance  problem.  In  this  case  it  may  be  difficult 
to  express  the  radiation  condition  in  moving  coordinates  in 
a  tractable  form. 

A  computer  program  can  probably  never  be  made  as 


efficient  or  as  general  as  possible.  However,  we  have  made 
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considerable  effort  to  attain  such  a  goal.  As  mentioned 

in  Chapter  IV,  our  computer  program  solves  the  coupled 

equations  of  cosine-mode  potential  and  sine-mode  potential 

simultaneously.  Accordingly,  the  band  width  of  the  final 

coefficient  matrix  is  double  the  band  width  of  the  either 

. —  „  .  — . -  — 

potential  of  either  mode.  It  would  be  worthwhile  to  try 
to  decouple  them  and  to  solve  them  separately.  This  might 
result  in  a  saving  of  computer  time. 

It  is  possible  to  extend  this  method  to  a  general  three- 
dimensional  boundary  geometry  without  much  difficulty. 

However,  the  band  vidth  of  the  coefficient  matrix  will  then 
be  very  large.  Hence  the  computing  time  will  also  be  very 
large.  If  we  modify  the  present  computer  program  slightly, 
we  can  solve  diffraction  problems  with  an  axi-symmetric  body 
in  three  dimensions  and  a  plane  incoming  wave.  Other  modes 
of  motion  than  heave  can  also  be  treated  for  the  axi-symmetric 
body. 
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APPENDIX  A 


The  potential  £(x,  y,  a,  b,  t)  due  to  a  source  with 
pulsating  strength  c®s  ft  at  (a,b)  is  given  in  (2-14),  i.e., 


Ianb,t)  *4  J  — K=v - ^J4**4* 

}. 


Sen  G'fc 

where  c  =  a  +  ib,  c  =  a  -  ib,  e  =  x  +  iy. 

Let  us  consider  a  contour  integral 
the  complex  K-  plane  (»8  shown  below) 
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From  the  cauchy  integral  theorem,  we  have 
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From  (A-2)  and  (A-3)  we  obtain: 
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In  order  to  make  the  computation  of  I  and  I  easy,  we 
shall  try  to  make  — ♦  0  as  l^-*0  and  shall  choose  the  path 

of  integration  in  I3  along  from  r  =  0  to  r  =  oo  so  that 

the  integrand  becomes  non-oncillatory  in  r.  From  this  we  obtain 
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(A -11) 


It  is  of  interest  to  note  two  simple  cases,  i.e.  and 

/«•  y 

(^-("^■3)  “J  <e2  =  °’  along  x_axis  x  *0)  (A-12) 

c°* 

SfLjfh)  “733^ -  (e2  =-f*  alon*  y  axls  y<0)  (A"13) 

From  (A-12),  (A-13)  we  can  see  the  integral  decrease  as  (V£,J  1 
along  the  x-axis  and  along  the  y-axis  for  large 

6>&) 

In  order  to  compute  (A-II),  one  may  use  the  Gauss-Laguerre 
formula  whenever  the  error  bound  for  the  numerical  integration 
is  permissible.  For  completeness  this  well-known  formula 
(Abramowitz  and  Stegun,  1967)  is  stated  here: 

j  j*>  e'XJx  =  |  «C  i<*£J  +  K*. ,  (A.14) 

Jb 

x^  is  the  i-th  zero  of  the  Laguerre  Polynomial  JL-*. (X) 
r  __  frl  / ) 1  ** _ 


where 
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APPENDIX  B 


We  shall  give  a  very  brief  description  of  general  quad¬ 
rilateral  elements.  A  reader  can  find  much  more  detail  in 
Zienkiewicz  (1971,  pp.  103-170)  and  Ergatoudis,  Irons  and 
Zienkiewicz  (1968).  The  procedure  for  the  numerical  inte¬ 
grations  in  (4-12),  (4-13)  and  (4-14)  will  also  be  described 
briefly. 

Let  ui  consider  a  four-node  quadrilateral  element  in  the 
physical  plane.  We  introduce  a  new  coordinate  system  in  such  a 
way  that  the  element  in  the  xy  plane  maps,  one-to-one,  into  a 
square  in  the  new  coordinates  (  ^  _  plane)  shown  in  Fig,  A-l, 

Then  we  can  express  a  function  in  terms  of  a  set  of  the 
interpolation  functions  [N]  ,  which  are  functions  of  ^  and  ^  . 

+  (B-l) 


j>  &  t  tJz  f. 


Fig.  A-l 

Similarly  we  can  write 


X  **  Mi  Z,  tx  +  ll» 
2  >» 
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where 


=  a~x>o-+t)/4  , 

7 

=  at  $)('-?)  A  ; 


(B-3) 


[N]  contains  the  interpolation  functions  In  the  new  coordinates. 
Hence  this  may  be  different  from  [N]  in  the  physical  plane. 

When  we  have  a  general  quadrilateral  element,  it  is  more 
convenient  to  perform  the  integrals  in  (4-12),  (4-13)  and  (4-14) 
in  the  new  coordinate  system.  From  the  coordinate  transformation, 


we  have 
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or  in  matrix  notation 


in  which  [J]  is  the  Jacobian  matrix, 
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The  integration  with  respect  to  x  and  y  in  (4-12)  can  simply 
be  changed  to  integration  with  respect  to  ^ and  y  ,  with  a 
simplification  of  the  limits  of  integration,  which  now  are  simply 
from  -1  to  l  in  both  variables,  and  with  the  change 


<**<3  =  \j\  j 
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where  jjJ  is  the  determinant  of  [J],  Similarly  the  boundary 
integrals  defined  in  (4-13)  and  (4-14)  can  also  be  carried 
out  in  the  new  coordinates  with  a  change  of  variables, 

or  _  _ 

-l  5  ~  JTxty)  ^ 

with  limits  from  -1  to  1. 

All  the  integrations  are  performed  numerically  using  the 
well-known  Gaussian  quadrature  formula  (Abramowltz  and  Stegun, 
pp.  916-917). 

In  our  computer  program,  the  element  shapes  that  we  used 
were  four-node  and  eight-node  quadrilateral  elements.  The 
procedure  for  computing  Kij,  Hij  and  bi  in  (4-12),  (4-1$  and 
(4-14)  for  the  case  of  an  eight-node  quadrilateral  element  is 
very  similar  to  that  for  a  four-node  element.  Therefore  this 
will  not  be  given  here. 

In  writing  the  computer  program,  a  computer  program  made 
for  structural  problems  (Wilson,  1970)  has  been  helpful. 
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Fig.  26b  -  Phase  Angle  and  Maximum  Wave  Amplitude 
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Fig.  27d  -  Phase  Angle  and  Maximum  Wave  Amplitude  of  the  Tot^l  Waves 
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Fig.  32a  -  Diffracted  Waves  Due  to  Two  Vertical  Flat  Plates 
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Fig.  33a  -  Diffracted  Waves  on  the  Free  Surface  (Two-Fluid  Problem) 
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